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Abstract: We present methods to compute higher orders in the threshold expansion 
for the one-loop production of a Higgs boson in association with two partons at hadron 
colliders. This process contributes to the N 3 LO Higgs production cross section beyond the 
soft-virtual approximation. We use reverse unitarity to expand the phase-space integrals in 
the small kinematic parameters and to reduce the coefficients of the expansion to a small set 
of master integrals. We describe two methods for the calculation of the master integrals. 
The first was introduced for the calculation of the soft triple-real radiation relevant to 
N 3 LO Higgs production. The second uses a particular factorization of the three body 
phase-space measure and the knowledge of the scaling properties of the integral itself. Our 
result is presented as a Laurent expansion in the dimensional regulator, although some of 
the master integrals are computed to all orders in this parameter. 
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1. Introduction 


The discovery of the Higgs boson [1] opens the path to the exploration of the sector respon¬ 
sible for the breaking of the electroweak symmetry and makes the Standard Model a fully 
predictive theory, because all couplings are now uniquely fixed. Precision measurements of 
the properties of the Higgs boson therefore represent a crucial test of the Standard Model, 
and any deviation, however small, will open a window to new physics scenarios. 

The main production mechanism of a Higgs boson at hadron colliders is gluon fusion. 
Unfortunately, gluon-fusion suffers from large perturbative instabilities [2, 3, 4, 5], and the 
uncertainty on the next-to-next-to leading order (NNLO) prediction still remains of the 
order of 10%. In order to fully exploit the potential of the LHC, it is therefore crucial to 
improve on the NNLO result by going to the next order in the perturbative expansion. 

Recently, the next-to-next-to-next to leading order (N 3 LO) QCD correction to Higgs 
boson production has been computed as an expansion around threshold, and it was shown 
that the remaining QCD scale uncertainty is reduced to a small 2% [6]. Most of the 
ingredients and subprocesses that go into the computation of ref. [6] have been published 
separately over the last few years: the three-loop corrections to Higgs production in gluon 
fusion, as well as the corrections from the emission of an additional parton at one or two 
loops, have been computed in full generality in refs. [7, 8, 9, 10, 11]. In order to obtain 
a finite result, appropriate ultra-violet and infrared counterterms need to be included [12, 
13, 14], The contributions from the emission of two partons at one loop and three partons 
at tree-level, however, have not been computed for arbitrary kinematics so far, but they 
are only known as an expansion around threshold. In particular, in refs. [15, 16, 17, 
18, 19] the first two terms in the threshold expansion of the triple-real and double-real- 
virtual corrections have been computed by reducing the corresponding phase-space and 
loop integrals to a small set of master integrals, all of which belong to a special class 
of integrals that we call soft integrals. These soft integrals are not only important when 
computing the first two terms in the threshold expansion of the cross section, but they also 
contribute to the full result for the cross section as boundary conditions to the differential 
equations satisfied by the master integrals in full kinematics. Despite their importance, the 
results for the individual soft integrals that contribute to the soft-virtual and next-to-soft 
corrections at N 3 LO of refs. [16, 19] have never been published explicitly 1 . 

The purpose of this paper is threefold: first, we close the aforementioned gap and we 
present all the technical details that went into the computation of the N 3 LO cross section 
in the soft-virtual and next-to-soft approximations of refs. [16, 19]. We discuss in detail 
methods to perform the threshold expansion of one-loop matrix elements for the production 
of a heavy colorless state in association with two partons to any desired order, at least in 
principle. Second, we present techniques to compute the resulting soft master integrals 
analytically as a Laurent expansion in the dimensional regulator with coefficients that are 
polynomials in multiple zeta values. Finally, we give explicitly all the soft master integrals 

1 Some of the soft integrals contributing to the soft-virtual corrections to the cross section at N J LO have 
been published in ref. [17] after the computation of ref. [16] had been completed. 
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contributing to the results of refs. [16, 19], as well as several new soft master integrals 
which only contribute at higher orders in the threshold expansion [6]. 

In a nutshell, our approach to threshold expansion can be described as follows: we 
start by expanding both the phase-space measure and the interference amplitudes into a 
power series in some small parameter quantifying the deviation from threshold. At every 
order in the expansion, the phase-space integrals are mapped onto (cut) Feynman integrals 
by virtue of reverse unitarity [3, 20], which identifies on-shell phase-space constraints with 
cuts of Feynman propagators. The resulting cut integrals are then reduced to a small 
set of soft master integrals using integration-by-parts (IBP) techniques [21]. Similarly, 
the expansion of the one-loop matrix elements in the integrand is performed using the 
strategy of regions [22], which allows to exchange the threshold expansion and the loop 
integration, provided that contributions from all ‘regions’ in loop momentum space that 
lead to singularities in the soft limit are included. Similar to the findings of ref. [8], we 
observe that in the present case the relevant regions correspond to regions where the loop 
momentum can be either ‘hard’, ‘soft’ or ‘collinear’ to one of the initial-state momenta. 
While the IBP reduction of the hard and soft regions can be dealt with using standard 
techniques, the collinear regions require some special attention, because the resulting in¬ 
tegrals contain propagators which are non-linear in the Lorentz invariants, and they are 
hence not amenable to standard techniques. We find that, due to the analytic structure 
of these propagators, it is always possible to recast the collinear-type integrals into a more 
standard form using partial fractioning, and we observe that the ensuing master integrals 
are always soft integrals, independently of the region. 

In a second part of the paper we discuss techniques to evaluate the soft master integrals 
as Laurent expansions in the dimensional regulator, and we present two complementary 
approaches to achieve this. The first technique is reminiscent of the method of ref. [15] 
and allows one to obtain a multifold Mellin-Barnes representation with poles at integer 
values for each soft master integral. The second method uses a particular factorisation of 
the three-body phase-space measure, similar to the method presented in ref. [23] for purely 
real emissions. This method allows to arrive at a simple parametric integral representation, 
whose simplicity can be understood via the scale invariance properties of soft integrals. 

The paper is organised as follows. In Section 2 we describe the setup of the problem 
and introduce some notation. In Section 3 we outline the methods used for the expansion 
around threshold and the technique of reverse unitarity. We discuss the reduction to 
master integrals and list the complete set of master integrals for each region in Section 4. 
In Section 5 we describe the methods used for their evaluation and we present analytical 
results, either as an expansion in the dimensional regulator or exactly in terms of generalised 
hypergeonretric functions. We draw our conclusions in Section 6. 


2. Double-real-virtual corrections to Higgs production 

In this paper we consider the production of a Higgs boson (H) in association with two 
massless partons (k and l ) via the scattering of two massless partons in the initial state (i 
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and j ) 


i{qi)+j(Q2) ->• #(<?#) + %3) + J(<?4) • 


( 2 . 1 ) 


Here the partons i,j,k and l can be gluons or massless quarks of Np different flavours, 
which are not directly coupled to the Higgs boson. We work in the large top-mass limit, 
and we assume the gluons to couple directly to the Higgs boson via the effective operator 

C es = -\CHG%G^. (2.2) 


At N 3 LO all of these partonic processes enter the inclusive cross section in the large top- 
mass limit through the interference of the relevant tree and one-loop amplitudes. To be 
more concrete, let us write the contribution to the cross section as 


<7 


RRV 


( s,z,e ) 


•A/ 'ij;kl 
2s 


/m(l) 


(z,e), 


(2.3) 


where Mij-ki contains all averaging and symmetry factors and s = s \2 denotes the square of 
the partonic center of mass energy. The mass of the Higgs boson will be denoted by rripr , 
and for later convenience we introduce the following notations, 


s h—ik — 




ik ’ 


Qh-ik = T hQii + ■ ■ ■ + TikQik 


where 

_ / +1 if i = 1,2, 

[-1 if i > 2, 

and 

Q 2 TYb TT _ 

QH = m H , 2 = -, Z = l-Z. 

s 

The dimensionless (coefficient) function C^\ klH is defined as 

Cij\klH( z ’ e ) = J 3 ^ 2Re (^ij\ H kl^ij\Hkl) ■ 

spins,colors 


(2.4) 


(2.5) 


( 2 . 6 ) 


(2.7) 


If we set the scale introduced by dimensional regularisation // equal to mu (a convention 
which we shall adopt throughout this paper), then the coefficient function only depends 
on z and the dimensional regulator e, related to the dimension of space-time D via e = 
Here A^\ Hkl denotes the L-loop amplitude for the the process ij —>• Hkl and the 
differential three-particle phase-space measure is defined as 

d$ 3 = {2ir) D 5 {D \ qi 234 ~ Qh) ^yD-i 5 +^ ~ II ( 2 ^d-i Mgi) * ( 2 - 8 ) 

where S + (q 2 — m 2 ) = 5(q 2 — m 2 ) @(g°), with 0 the Heaviside step function. In the rest 
of this paper we present methods to compute the coefficient function C^\^ klH (z) as an 
expansion around threshold. 
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3. Threshold expansion for real emissions at one loop 


Our goal is to obtain the expansion of the coefficient function C^\^ klH (z,e) close to z = 1, 
where all the emitted final-state partons are soft. In other words, we want to expand 
eq. (2.7) in the soft momenta q 3 and ( 74 . To be more concrete, let us define a set of rescaled 
momenta pi as 

_ / zpi, if i € {3,4} , 

Qi — { ,1 • ( 3 - 1 ) 

I pi , otherwise . 

After this rescaling both the phase-space measure and the matrix element in eq. (2.7) 
depend on the scaling parameter z. It will be convenient to introduce a notation for how 
a given quantity scales with z, e.g., we will write q\ ~ 1 and 1/3 ~ z to refer to the scaling 
behaviour defined by eq. (3.1). Our goal is to expand both the phase-space measure and the 
amplitudes in this scaling parameter. In the rest of this section we discuss the expansion 
of each of these objects in turn. 

Let us start by discussing the threshold expansion of the phase-space measure. Chang¬ 
ing variables to the rescaled momenta (3.1), it is obvious that the only quantity in eq. (2.8) 
that depends non-trivially on the scaling parameter is the on-shell condition for the Higgs 
boson, 


S + (p 2 H -m 2 H ) = 5 + ((p 12 - zpu) 2 - m 2 H ) =-5 + (s 12 - 2p 12 ■ p 3 4 +zs 34 ). (3.2) 

At leading order in the threshold expansion we can simply ignore the term linear in z 
appearing inside the <5 function. If we want to study the subleading terms in the ex¬ 
pansion, however, we need to take into account this term, i.e., we need to expand the 5 
function around threshold. This can be achieved using reverse unitarity [3, 20] to interpret 
phase-space integrals as Feynman integrals with cut-propagators, i.e., where some of the 
propagators have been replaced by on-shell (5-functions, 


5+(q 2 ) 



(3.3) 


The subscript c is a reminder that this propagator is cut. Cut propagators can be differ¬ 
entiated just like normal propagators, 


d / 1 \ n _ / 1 \ n+1 d F(x) 

d x\F(x)J c n \F(x)J c dx 

but satisfy the extra condition 


(q z ) m = { W 

0 , 


if n > m 
if m > n 


The soft expansion of the phase-space measure therefore reads, 


d $ 3 = $f’ n) z n , 

n=0 


(3.4) 


(3.5) 


(3.6) 
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with 


dd> 


(S,n) 

3 


2vr (-s 34 ) n 


\ l+n 4 a £) 

1 \ TT & Pi 

S12 - 2^12 ' f>34 / c ^ (27r) Z) - 1 



(3.7) 


Note that this procedure necessarily introduces cut propagators raised to higher powers in 
the subleading terms in the expansion. In Section 4 we will see that integrals involving 
these additional powers of cut propagators can always be reduced to the case n = 0. As a 
result, all phase-space integrations will be performed against the soft phase-space measure 
defined by 

d4>f = d$f’ 0) = 2 tt<5 + (si2 - 2 p 12 • p u ) h(Pi) ■ (3-8) 


Let us now turn to the threshold expansion of the integrand. This requires expanding 
in z both the tree-level and one-loop amplitudes appearing inside the matrix element in 
eq. (2.7). The expansion of the tree-level amplitude can easily be obtained by introducing 
the rescaled momenta (3.1) and expanding the resulting rational function in z. 

The expansion of the one-loop amplitude, however, is more subtle. In order to under¬ 
stand why, let us mention already now that the coefficient function is not meromorphic at 
z = 1, but it can nevertheless be written in the form 


y —*( 1 ) 


(*.<0 


= z- 1 ^' 


C ij^llH( z ’ € ) + z 1 beC ij^klH( Z ’ e ) + 2 1 beC ij^’klH( Z i e ) 


1 —56^(1,c) 


1—6e/'-f(l>s) 


(3.9) 


where the functions C^f^ klH {z) for r G {h,c,s} are holomorphic at z = 1 and therefore 
admit a Taylor expansion around this point, 


r^ 1 ’ r ) r z 


e = 


E c , 


(1 ,r,k) 

ij^tklH 


,e z 


k =0 


(3.10) 


A naive expansion of the loop integrand in z fails to reproduce this analytic structure, 
because the expansion of the loop integrand gives rise to a holomorphic function, and the 
phase-space measure can only account for the term proportional to z~ 4e (see eq. (3.6)). In 
other words, the threshold expansion does not commute with the loop integration. 

Indeed, it is well-known that for some kinematic limits the loop integration does not 
commute with expansions of the loop integral, and hence a naive Taylor expansion of the 
loop integrand is bound to fail. The strategy of regions instead allows one to recover the 
correct expansion of the integral [22], In this approach one sums up the expansions around 
all loop momentum regions which can be identified to lead to singularities in the limit z —» 1. 
In ref. [8] it was argued that in the present case these different ‘regions’ can be classified 
as ‘hard’ (h), ‘soft’ (s) and ‘collinear’ (c), depending on whether the loop momentum is 
hard, soft or collinear to one of the hard partons. The three different terms in eq. (3.9) 
are then associated with these three regions. The collinear region is further decomposed 
into two different collinear regions, ‘collinear-1’ (c 4 ) and ‘collinear-2’ ( 02 ), where the loop 
momentum is collinear to p\ or p 2 - Each region is associated with a specific scaling of 
the (components of the) loop momentum k, similar to the rescaling (3.1) used to define 
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the threshold expansion of the phase-space measure. In order to define these scalings, we 
parametrise the momentum flowing through a suitable propagator by 



k^ = a pi + /3 P 2 + kj_ , 

(3.11) 

where 


k ■ P 2 a k-pi 

a = -, p =- 

(3.12) 

and the 

Pi -P2 Pi' P2 

measure can then be written in this parametrisation as 


& D k = da d/3 d D ” 2 A:_L . 

(3.13) 


The different regions are then associated with the following scalings in z, 

(h ): k 1 , such that d D k 
(ci): a ~ z, k± ~ y/I, such that d D k ~ fa", 

(02): j3 ~ 2, k± ~ VI, such that d^A; ~ Fa, 

(s): A: ~ z, such that d D A: ~ z D . 

The overall scaling of a given region in eq. (3.9) can then be understood by combining the 
scalings of phase-space measure, eq. (3.6), and loop momentum measure in the different 
regions. 

At this stage we have to address a critical point in this procedure. The scalings we 
have just defined are not invariant under the shift symmetry of the loop integral, i.e., 
the rescaling can only be applied once an adequate shift of the loop momentum has been 
identified. We solve this issue in the following way: first, we note that, up to swaps of the 
external momenta ^ 3 ,^ 4 , all the loop integrations appearing inside the one-loop amplitude 
can be mapped to one of the following four pentagon topologies, 

Pent mj ... >n 5 (pi,-p 3 ,-p 4 ,P2), 

Pent m (pi,p 2 , -PA, ~P'i), 

Pent ni! ... )n5 (pi4,p4,p 2 4, -P3) , 

Pent ni! ... >n 5 (-p 3 ,pi,p 2 ,-P 4 ), 

with 

Pent ni) ... )n6 (gi, 92,93,94) 

_d D _k _ 1 _ 

iir D / 2 [(k - 91 - 9 2 ) 2 ] ni [(A; — 92 ) 2 ] n2 [A: 2 ]” 3 [(A; + 93) 2 ] n4 [{k + 93 + 94) 2 ] n5 ’ 

Using the loop momentum scalings defined in the previous paragraph we have checked 
explicitly that the first few terms of the z expansions, obtained by summing up the expan¬ 
sions of non-vanishing hard, collinear and soft regions of various pentagons, boxes, triangles 
and bubbles defined through eq. (3.14), reproduce the same expansions as one can obtain 



(3.14) 


(3.15) 
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from an expansion in Feynman parameter space by the methods of ref. [24]. Combined with 
the expansion of the phase-space measure discussed at the beginning of this section, we 
have therefore obtained a machinery to compute the threshold expansion of the coefficient 
functions c\]\ klH (z,e). 

Let us conclude this section by making a technical comment about the expansion in 
the collinear regions. The procedure outlined above will effectively lead to an expansion 
in sfM in the collinear region, because the transverse components of the loop momentum 
scale like y/2. hi other words, there seems to be a contradiction to the statement that the 
coefficient function in the collinear region, C^^ klH (z,e), is holomorphic at 2 = 1. In the 
following we show that this coefficient function in the collinear region is indeed holomorphic 
at z = 1, and therefore all terms proportional to powers z n / 2 , with n odd, must vanish. In 
order to prove this statement, consider an integral of the form, 


r d D k A^Af ...Ay 

J i7T D / 2 F(k,pi,p 2 ) 


(3.16) 


where F is a Lorentz invariant function which only depends on the loop momentum k 
through the scalar products A 2 , k-p\ and k-p 2 . Using the scaling properties of the collinear 
regions and the relations of eq. (3.12), it can be seen that all integrals which appear in the 
z expansions of the collinear regions of the four pentagon topologies defined in eq. (3.14) 
indeed fall into this category. As the only source of in the expansion are the transverse 
components of the loop momentum Aj_, all terms in the expansion proportional to z n / 2 
with n odd must be proportional to an integral of the type (3.16). A sufficient condition 
to prove that C-j^ klH (z, e) is holomorphic at z = 1 is therefore that all integrals of the 
type (3.16) vanish for odd values of n. In Appendix A we prove that this is indeed the 
case. More precisely, we prove the following result, 


r d D k A^Af ... 

J in 0 / 2 F(k,pi,p 2 ) 

UlH2---Hn _ _ _ 

= < x nri?(£-4+2 i) 

0 


d D k (Aj_) n / 2 
m D ! 2 F(k,pi,p 2 ) 


if n even, 
if n odd, 


(3.17) 


where 


ffT = r 


PlP2 + p\p2 


Pi 'P 2 

is the metric on the space transverse to the vectors p\ and p 2 , and 


UlH 2 -Hn 

y± 


1 

2”/ 2 (n/2)! 


gi^( 1 )^( 2 ) 


9± 


crCSn 


Mcr(n— 1) Mcr(n) 

■■■ 9 _ L 


(3.18) 


(3.19) 


We note that eqs. (3.17) and (3.11) can be used efficiently to reduce any one-loop integral, 
which appears in the collinear region, containing arbitrary powers of the Lorentz invariants 
A-p 3 and A-jq in the numerator, to integrals containing only the Lorentz invariants k-pi, k-p 2 
and A 2 in numerator and denominator. 
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4. IBP reduction and master integrals 


In the previous section we argued that the threshold expansion of the coefficient function 
cjj\kiH receives contributions from hard, soft and collinear regions, see eq. (3.9). The 
integrals appearing in the expansion are, however, not independent, but can be reduced to 
a small set of master integrals using integration-by-parts (IBP) identities [21]. While IBP 
identities have been introduced for Feynman integrals rather than phase-space integrals, 
we can use reverse unitarity to interpret phase-space integrals as Feynman integrals with 
cut-propagators, to which IBP identities are known to apply. Here we go one step further, 
and we apply IBP identities combined with reverse-unitarity to the individual hard, soft 
and collinear regions. As a result, we obtain for each region a small set of hard, soft or 
collinear master integrals in terms of which the coefficients in the threshold expansion can 
be expressed. Independently of the region, we will see that the master integrals will fall 
into the class of so-called soft integrals, i.e., integrals with respect to the soft phase-space 
measure (3.8) of a function that is independently homogeneous under a rescaling of the 
initial momenta p\ and p 2 , as well as under a simultaneous rescaling of all the final-state 
soft momenta. In the remainder of this section we review the IBP reduction in each region, 
and we present the analytic results for the master integrals in each region. Details about 
the computation of the master integrals will be given in Section 5. 

4.1 The hard region 

We start by discussing the master integrals coming from the hard region. Since the loop 
momentum is hard, we can immediately expand in the soft final state momenta under the 
integral sign, and the IBP reduction for the combined loop and phase-space integral follows 
the same lines as for the purely real soft emission discussed in ref. [15]. We find that the 
whole contribution from the hard region can be expressed in terms of only two master 
integrals, 


M? 





s 


d$f Bub(si 2 ), 



Bub(s 12 ). 

Sl3 s 24 S34 


The double line denotes the Higgs boson, and the dashed line represents the phase-space 
cut. All other internal uncut lines are scalar propagators. Moreover, Bub(p 2 ) denotes the 
usual one-loop bubble integral, 



d D k 1 

iir D / 2 k 2 (k + p ) 2 


cy 

e(l-2e) 


(-P 2 - i 0 )~ £ . 


(4.1) 







where we defined the usual loop factor 


r(l + e)r(l -e ) 2 
r(l - 2e) 


(4.2) 


The computation of the master integrals will be discussed in detail in Section 5. Note 


that Bub(s) = Bub(m^) + O(z), and so the hard master integrals can be written as a 
bubble integral Bub(m^) multiplying one of the two soft master integrals appearing in the 


double-real soft corrections at NNLO. 

4.2 The soft region 

Next, we discuss the contribution from the soft region, where k ~ z. It is straightforward 
to extend the methods of ref. [15] to the soft region at one loop. Indeed, since the loop 
momentum scales in the same way as the final-state soft momenta, IBP reduction proceeds 
in the same way as for purely real soft emissions at tree-level, except that the loop momen¬ 
tum is not constrained to be on shell. At the end of this procedure, we find the following 
set of soft master integrals, 






in D / 2 (~ 2 kpi)( 2 kp 2 - 2p 2 P3 - ^P 2 Pi){k 2 - 2kp 3 ) ’ 
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1 



1 
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d D k 1 

d<ho -;- 

3 m 0 / 2 ( k 2 — 2kpz — 2kp4 + 2p^pi)k 2 
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M f 



_i_, 

J 3 m D ! 2 ( k 2 + 2 kp 3 )(k 2 - 2kp 4 ){2kp 2 - 2p 2 p 4 ){-2pip 3 - 2kp\) ’ 



/ d$ f^L_I_, 

J 3 in D / 2 ( -2p 2 p 3 )k 2 (2kp 2 - 2p 2 p 4 )(-2p 3 p 3 - 2kpi)(k 2 + 2 kp 3 ) 



/d*!-^._I___. 

J 3 zvr 15 / 2 ( -2p 1 p 3 )(-2p 2 p 4 )(k 2 - 2kp 4 )(k 2 + 2 kp 3 ) 


Ml 



f ,*s d ° k 1 

/ d$Q- 
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/«_1_ 
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1 

(2P3P4 ) ( - 2pip 3 ) (■ -2p 2 p 4 ) 


M s n = 
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_I_, 

J 3 m D ! 2 ( k 2 + 2kp 3 )k 2 (2kp 2 - 2p 2 p A )(k 2 - 2kp A )(-2pip 4 ) ’ 


Mf 2 
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M s l4 = 

1 -*■ 

= td*§*°* _i_. 

J 3 iTT D/2 k 2 (-2kp 1 )(2kp 2 - 2p 2 p 4 )(k 2 - 2kp 4 )(2p 3 p 4 )(2kp 2 - 2p 2 p 3 - 2p 2 p 4 ) 



The number of propagators of the soft master integrals depending on the loop momenta 
ranges from two to five. Since we are working in the soft region, the loop integrals are not 
the complete momentum space integrals, but they correspond to the eikonal approximation. 
Note that the eikonal approximation of the loop integrals depends on the specific orienta¬ 
tion with which it enters the soft phase-space integrals, or equivalently, which invariants 
become soft. We have explicitly evaluated all loop integrals in the eikonal approxima¬ 
tion by determining the scaling of the Feynman parameters in the limit using the package 
asy.m [24], We find that in all cases, except for one pentagon integral, the remaining 
parametric integrations can be perfomed in closed form to all orders in the dimensional 
regulator e. For the remaining pentagon we are able to obtain Mellin-Barnes integral repre¬ 
sentation valid to all orders in e. In the following we give a brief summary of the results for 
the soft virtual integrals. The bubble type integrals in A4 3 and Mf can just be computed 
as in eq. (4.1), and will not be discussed any further. Note that we present the results in 
the Euclidean region where all invariants are negative, and the analytic continuation to the 
scattering region is given by 

—m 2 H — iO^-e~ in m 2 H , —si 2 — i0^e~ m si 2 , —s 34 — iO^-e~ m s 34 . (4.3) 

Note also that the loop integrals appearing inside the master integrals have the correct 
homogeneity properties to turn the master integrals into soft integrals. The computation 
of the master integrals will be discussed in Section 5. 


Soft Triangle Integrals. The soft master integrals A4f and A4f contain a virtual in¬ 
tegral with three propagators. 

-P1234 



r(l-e)r(l + e) 2 


(—S 13 ) e ( — s 24 ) e (-si 2 ) 


-l+e 


The integral corresponds to a completely off-shell triangle where two of the virtualities of 
the external momenta (s 24 = p% 4 and S 13 = p\ 3 ) scale as z. 
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Soft Box Integrals. We observe that there are four configurations of loop integrals with 
four propagators. Two of them arise from box master integrals with one off-shell leg and 
two from configurations with two off-shell legs. 

1 . The loop integral appearing in the definition of Adf is given by 


P3 


PI 



P2 


dPk 


m D ! 2 k 2 {k — P3) 2 (—2kpi)(—2p2P3 + 2kp2) 


~P 123 


2(- 5 i2) e (-si3)- £ " 1 (-S23r £ - 1 r(l - e) 3 T(e + 1)- 
e 2 T(l - 2e) 


2 . The soft master integrals A4f, Ad fj and Adf 2 contain the loop integral 




d D k _1_ 

m D ! 2 k 2 (k — P3) 2 (k — P3 — p/i) 2 (—2kp2) 


2 (-s 3 4 ) _e ~ 1 r(l - e) 2 T(e + 1) „ (, 1 , *24 \ 

e 2 (s23 + s 2 4)r(l-2e) 2 1 V ’ ’ ’ ^23 + W ’ 

where 2 T 1 denotes Gauss’ hypergeometric function, 


OO 

2 -F 1 (a, b; c;z) = ^ 

n=0 


{fl)n (^)n 

(c)n 


Z 

n\ 


(4.4) 


3. In A4* and we find a limit of a two-mass-easy box integral. 
P34 Pi 

r d D k 1 



in D / 2 {k + p 3 ) 2 (k - Pi) 2 { 2 kp 2 — 2p2PA){—2piP3 — 2kp\) 


P1234 


= ( '-’ 34 ) 


2T(1 — e) 2 T(e + 1) 


(e + l) 2 (si 3 + * 14)(*23 + s 2 4)r(l — 2e) 

S12S34 


x 2 F 1 I 1,1; e + 2; 1 — 


(*13 + *14)(*23 + s 2i) ) 


- 13 - 














4. The four propagator loop integral appearing in master integral Atf is a limit of a 
two-mass-hard box integral. 


P3 


P24 




d D k 1 

m D ! 2 k 2 (k + p 3 ) 2 (-2pip 3 — 2kpi)(2kp 2 - 2p 2 p 4 ) 


= (-si 2 ) e (-si 3 ) 


« 1(s24) « 1 r(i- i )r( 1± i) !2Fi 


1, e + 1; 1 — e; — 


£23 

5 2 4 


The loop integral in A4g is related to the integral above via the permutation of 
external momenta (p\ —> p 2 ,p 2 —>• Pi,P a —> p 3 ,p 3 —>• p 4 ). 


Soft Pentagon Integrals. We also need the soft limits of two configurations of the 
one-loop pentagon with one off-shell leg. 


1. Unlike the other one-loop soft integrals we require, the pentagon integral appearing 
in At f 3 does not permit a closed expression in terms of simple 2 T\ hypergeometric 
functions. Nonetheless, we are able to derive a compact Mellin-Barnes representation. 


—P 1234 



/ 


d D k _1_ 

m D ! 2 k 2 (k + p 3 ) 2 (k - p.i) 2 {2kp 2 — 2p 2j p i )(—2pip 3 - 2kp \) 


2e(2e + 1)(—fi 34 )~ 6 ~ 2 f +i °° dz x dz 2 d^ 3 

si 2 T(1 - 2e) J_ ioo (2vri) 3 


£13 + 514 V 3 f 5 24 V 2 /£l 3 (£ 23 + 5 2 4)y i 
513 / \ 5 23 + S 24 / V 5i 2 S 3 4 / 


xr(zi + i)r(-z 2 )r(-z 3 )r(z 2 - zi)v{z 3 - £i)r(2i + e + 2) 


xf(-z 2 - e - 1) T(«i + z 2 - z 3 + l)T(-zi + z 3 - e - 1). 
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2. The pentagon integral required for -Mf 4 is given by 


—P1234 



d D k 1 

iir D / 2 fc 2 (-2pip 4 - 2kp 1 )(2kp 2 - 2p 2 p 3 )(2kp 2 )(k + p 4 ) 2 


= s 42 (-si 4 ) e 1 (-«23) e 2 

+S 42 (-Sl4) _e_1 ( - S24) _e_1 


r(l-e)T(e + l) 2 , 

- 7 - 21 

2T(1 -e) 3 T(e + l) 2 
S23 e 2 T(l - 2e) 


1, e + 1; 1 — e; — 


£24 

•S23 


4.3 The collinear region 

In this section we discuss the reduction to master integrals in the collinear region. Unlike 
the soft and the hard regions discussed in the previous sections, the phase-space integrals 
in the collinear region cannot be reduced to master integrals using standard techniques. 
This is due to the appearance of propagators which are non-linear in the combined phase 
space and loop Mandelstam invariants. Effectively, the IBP relations do not close on certain 
topologies, and we are not aware of a method capable of dealing with IBP identities relating 
integral across different topologies. However, as we already pointed out in Section 3, 
all the one-loop integrals appearing in the collinear region have the property that their 
denominators only depend on the loop momentum k through the Lorentz invariants k 2 , k-p\ 
and k-p 2 . This brings as a consequence that there exist relations among the denominators, 
which can be used to find partial fraction identities to reduce all pentagons and boxes 
to a number of different triangle topologies. After integration over the loop momentum 
these triangles can be identified as linear combinations of only two types of bubbles. As a 
consequence, we can write the contribution from the collinear region in the form 




(*) 



Adsy}, e, 2 )Bub(si 3 ) + B({sij}, e, z)Bub(si 3 + s i4 ) , 


(4.5) 


where A and B are expanded as power series in z, whose coefficients are rational functions 
of the Mandelstam invariants and the dimensional regulator e. We dropped the real part 
of the loop integrals, because they are real for si 3 , si 4 < 0. We emphasise that eq. (4.5) 
is true to all orders in the expansion parameter z and can be proved simply by inspecting 
the propagators in the collinear region. 

Next, we express bubble integrals in terms of tadpoles via the identity 

Bl ,b( S ) = -cr , _ 2 e r (1 + e) a (4 ' 6) 

Since the single denominator in eq. (4.6) is clearly linear in all Lorentz invariants, this 
representation allows for a straightforward and efficient application of the IBP reduction 
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technique to the combined soft phase space and (now trivial) one loop integral. We find 
that in the collinear region we can reduce all integrals to just four soft master integrals, 


M? = 


M% = 


M<i = 


M c a = 

where the dotted lines represent numerator factors. We stress that the master integrals 
appearing in the collinear region are again soft integrals, and hence they can be computed 
using the same techniques as the master integrals in the soft and hard regions. This will 
be discussed in Section 5. 

5. Evaluation of soft one-loop integrals 

In the previous section we have seen that, independently of the region they originate from, 
the master integrals appearing in the threshold expansion are always soft integrals, i.e., 
they all take the form 

M\ = f d$f F t I (p 1 ,p 2 ,P3,P4-,e), (5.1) 

where d^g denotes the soft phase-space measure of eq. (3.8), and the integrand is homoge¬ 
neous with respect to a rescaling of the soft momenta as well as of the initial-state momenta 
Pi and p 2 , 

Fi(Xpi,P2,P3,P4 ! ;e) = \ au F l I (pi,p 2 ,P3,P4]e ), 

Fl(pi,Xp 2 ,P3,P4-,e) = F l I (p ll p 2 ,P3,P4-,e) , (5.2) 

f!(pi,P2, Ap 3 , Ap 4 ;e) = A 7i/ F, 1 (p 1 ,p 2 ,p 3 ,P4;e), 
for some an, fin, 7 * 7 . 



= / d$f (-S 13 - s u ) Bub(si 3 + S 14 ) 



= / d$ 3 (-si3)Bub(si 3 ), 




= / d$ 


S (~ s 13 ~ Sl4) 
Sl3'S24S34 


Bub(si 3 + S 14 ) 


(4.7) 



d $f Bub(s 13 ), 
S24S34 
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In this section we present methods to evaluate soft integrals analytically. We discuss 
two methods: the first is based on the derivation of a Mellin-Barnes representation for 
soft integrals, while the second exploits a particular phase-space factorisation as well as 
the homogeneity of the soft integrals to derive a parametric integral representation for 
soft integrals. For an additional method, using differential equation techniques, we refer 
to ref. [23]. We note that although we concentrate exclusively on the case of double-real 
emissions at one loop, all these techniques can in principle be generalised to arbitrary 
numbers of loops and legs in an obvious way. 


5.1 Mellin-Barnes representations for soft integrals 


In this section we describe a general method to write a given soft integral as a (possibly 
multi-fold) Mellin-Barnes (MB) integral with poles at integer values of the integration 
variables. In ref. [15] such a procedure was introduced for soft integrals of purely real 
emissions at tree-level. In this section we briefly review the procedure of ref. [15], and in 
the end we argue that it can easily be extended beyond tree-level. As the procedure is 
essentially identical to the purely real case, we will be brief and refer to ref. [15] for details. 

The procedure of ref. [15] starts from the observation that at tree level one may assume 
without loss of generality that the integrand F is a product of powers of two-particle 
Mandelstam invariants. Indeed, this can always be achieved by replacing every sum of 
two-particle invariants in the denominator by its MB representation, using the well-known 
formula 


1 

JaTW 


i 

m 



dz 


27TZ 


T r(-z)r(z + A) 


A z 

B z + X ’ 


(5.3) 


where the contour separates the poles at z = n from those at z = —A — n, n € N. 

Next, we parametrise the soft phase space using the energies and the angles of the soft 
momenta in the center-of-mass frame of the initial-state system. We write, with s \2 = 1, 


Pi = ^(M,0, •••)> P2 = ^(1,-1,0,...), pi = ^EiPi,i = 3,4, (5.4) 


where /% is the four-velocity in the direction of p t . The soft phase-space measure becomes, 


d$f = (2vr) 3 “ 2D 2 2 ~ 2D 5(1 -E 3 - E a ) E°~ 3 E°~ 3 d E 3 d E A dflf” 1 dfif' _1 , (5.5) 

where parametrises the solid angle of the soft momentum p{. Due to our assump¬ 

tion that the integrand is a product of powers of two-particle Mandelstam invariants, the 
integration over the energies is simply a Beta function. The remaining angular integrals 
can be written as MB integrals with poles at most at integer values of the integration 
variables [25, 26]. 

If we follow this procedure, every tree-level soft integral can be written as a multifold 
MB integral with poles at integer values of the integration variables. If the integrand 
contains loop integrals which evaluate to complicated special functions, this claim is no 
longer necessarily true. It does, however, stay true if the loop integral itself admits an MB 
representation of the same type. In our case, the integrand contains at most hypergeometric 
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(5.6) 


functions, which admit the MB representation 


2 -P 1 (a, b; c;x) 


T(c) f +ioo dz T(z + a)T(z + b) 

r(o)r( 6 ) J- ioo 2iri Z r(z + c) 


An exception is one pentagon integral, which we did not express in terms of simple hyper¬ 
geometric functions, but it still admits a multifold MB representation of a similar type, 
see eq. (4.5). We can therefore obtain an MB representation of the desired type for all soft 
integrals considered in this paper. 

The MB integrals obtained in this way can be evaluated using standard techniques. In 
some cases it is possible to perform all the MB integration in closed form without expanding 
in e, and one obtains generalised hypergeometric functions, 


OO 

pFq(a \,..., 0 ^; bi ,..., bq\ z) = ^ ' 

n=0 


n?=i(*i)» 



(5.7) 


Whenever we are not able to obtain a closed expression in terms of hyper geometric func¬ 
tions, we resolve the poles in e using standard techniques [27]. The result is a Laurent series 
in e whose coefficients are MB integrals whose contours are straight vertical line. In all 
cases these integrals can be evaluated numerically in a fast and efficient way. Alternatively, 
one can close the integration contours at infinity and sum up the residues of the poles of 
the Gamma functions in terms of nested harmonic sums [28] that evaluate to multiple zeta 
values. 


5.2 Soft integrals from phase-space factorisation 


We present in this section an alternative way to compute the soft integrals of Section 4, 
based on a factorisation of the phase space that separates the soft part of the phase space 
from the phase space for the emission of the Higgs boson. 

We start by writing the phase space for the production of a Higgs bosons in association 
with two massless partons as a convolution 

d $3 = J ^dT 2 (m^,^ 2 ;pi 2 )d4>2(0,0;Q), (5.8) 

where d4>2(mf, m|; q) denotes the phase-space measure for the decay of a heavy state with 
momentum q into two particles with masses mi and m 2 , 


d$ 2 (m?,m|;g) = (2tt) d 6^ D \q-q 1 -q 2 ) 6 + (qj-inj)8+(q, 


Note that for now we work with the full phase-space measure, and we expand 
end. A soft integral of the type (5.1) can then be written as 


ml). (5.9) 
in z at the 


M{ = j ^d$ 2 (m^,^ 2 ;p 12 ) J 7 !(pi,P2,Q) , (5.10) 

where we defined 

J 7 i(pi,P 2 ,Q) = J d$> 2 ( 0 , 0 ;Q)F i I (p l ,p 2 ,p 3 ,p 4 ;e ). (5.11) 
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Let us first concentrate on the computation of the integral J 7 /. The homogeneity of the 
integrand of the original soft integral, eq. (5.2), combined with the fact that the integration 
measure does not depend on p\ and P 2 , implies that T\ is homogeneous under a rescaling 
of any of its arguments, 

^/(Api,P2 ,Q) =\ aiI J 7 l I (pi,P2,Q), 

T\ (pi, Xp 2 , Q) = A fe T\ (pi,p2, Q ), (5.12) 

J 7 !(pi,P2,XQ) = A 7iJ J 7 l(pi,P2,Q)- 


Lorentz invariance then implies that the non-trivial functional dependence of can only 
be through the ‘cross ratio’ 


(gl • P-2)Q 2 
2 (pi • Q) (P2 • <3) ’ 


(5.13) 


Without loss of generality, we may write 


J 7 i{Pl,P2,Q) 


(Q2)(7U-<Xii-Pii)/2 i 

(pi ■ Q)~° u (p2 • <2) _/3i/ ^ ^ 


(5.14) 


We can also give a kinematical meaning to the cross ratio u: it is related to the angle 9 12 
between p\ and j <2 in the rest frame of Q, 


u = 


1 — cos 6 12 
2 


(5.15) 


Note that this cross ratio is precisely the argument of the hypergeometric function in the 
box integral appearing in A4§ and A4‘f 0 . 

This suggests that the most natural frame in which to parametrise the phase space in 
eq. (5.11) is the rest frame of Q. Writing pi = 4£)/% with fi t = (1 ,Uj), the phase-space 
measure in this frame becomes 


d$ 2 (0, 0; Q) = 2 1 ~ d (2? t) 2 ~ d (Q 2 )(-°- 4 )/ 2 dllf- 1 
and the invariants are 

(Pi ' Q)(P2 ■ Q ) 


(5.16) 


•534 = Q , S12 = 2 ■ 


Q 2 


(1 - COS0l 2 ) , 


S13 = ~{P1- Q) Pi' @3, «14 = ~{PI ■ Q) Pi ■ P3 , 

«23 = (p2 ■ Q) P 2 ■ P3 , «24 = -(P2 ‘ Q) P 2 ' P3 , 


(5.17) 


with Pi = (1,— Hi). Here we used the fact that in the rest frame of Q the final-state 
massless particles are back-to-back and their energy is \JQ 2 /2. The remaining integral 
over the solid angle can in all cases be performed using the formula [25, 26] 


dn 


D—l 


(Pi • Ps) m (Pj • PsY 


c\2—m—n—2e _1—e 
= Z 7T 


T(1 — m — e)T(l — n — e) 
T(1 — e)T(2 — m — n — 2e) 


2-^1 


^m, n; 1 — e; 1 



(5.18) 
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(5.19) 


The scalar product appearing in the right hand side can take the following values, 

Pi- P 2 = Pi ■ Pi = 1 - cos 012 = 2 u , 

Pi ■ P 2 = Pi ■ P2 = 1 + COS 012 = 2(1 - u). 

Note that in some cases the integrand may contain hypergeometric functions depending 
on angular variables. We can reduce the problem to the angular integral (5.18) by intro¬ 
ducing an MB representation for these hypergeometric functions. Eventually, we can in 
this way explicitly determine the function // (u) in eq. (5.14). We stress that although the 
computation was done in the rest frame of Q , the result is independent of the frame. 

Inserting eq. (5.14) into eq. (5.10), we are left with the computation of the remaining 
phase-space integral in eq. (5.10). We write Q = api+Pp 2 + Q±, withpi-Q_i_ = P 2 'Q± = 0 
and we work in the center-of-mass frame of the collision, and the on-shell conditions for Q 
and the Higgs boson enforce a + P = 1 + 0(z) and Q 2 ± = a( 1 — a) — p 2 . Note that this 
implies 0 < a < 1 and 0 < p 2 < a( 1 — a). In terms of this parametrisation the invariants 
are 


S 12 = 1, Q 2 = p 2 = a (1 — a) u , 2 pi • Q = 1 — a , 

and the phase-space measure can be written as 

,2 ,2_\ _ 1 tr,\-2+2e r ,/i „.wi „.M —e A A nD-2 


2 P2-Q = a, (5.20) 


d$ 2 {rriH,p 2 ;pi2) = ^ ( 27r ) 2+2 ' 


[a (1 — a) (1 — u)] e da d Qq 2 + 0(z ). 


(5.21) 


where ^q 2 parametrises the solid angle of Q in the center-of-mass frame. Putting every¬ 
thing together we see that the integrals over a and over the solid angle are trivial, and so 
we get, with 6 = (7 — a — P)/2, 

1 T (2 — e + a + <5)T(2 — e + /? + <5) rl 


M = 


1 2 a +P (4vr) 2 - e T(1 - e)T(4 - 2e + a + P + 26) 


duu 6 (l-u)- 6 fl{u). (5.22) 


We obtain in this way a representation for the soft integral as a simple integral over 
the function //. This last integral is usually easy to carry out. In many cases it can 
be performed using the Euler-integral representation of the generalised hyper geometric 
function, 


3 ^ 2 ( 01 , 02 , 03 ; h, b 2 ; z) 

r (0 2 ) 


dM «° 3 1 (1 — u) b2 0-3 1 2-01 (ai, 02 ; 61 ; uz). 


(5.23) 


r(o 3 )r(6 2 — 03) J o 

In those cases where the function fj is more complicated, we can either insert an MB 
representation for it, or alternatively perform the integral over u order by order in e in 
terms of harmonic polylogarithms [29]. 


5.3 Analytic results for the master integrals 

In this section we present the analytic results for all the master integrals introduced in 
Section 4. As only the real part of the master integrals enters the final result for the cross 
section, we only present result for the real part, which we normalise according to 

Re (Mi) = c r cos(Tre) 4>f (e) M\ , (5.24) 
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where cr is defined in eq. (4.2) and <3?f (e) denotes the soft phase-space volume, 



1 T(l-e) 2 

2 (4vr) 3 - 2e r(4 - 4ej ' 


(5.25) 


We have computed all the master integrals of Section 4 using the two different approaches 
described in the previous section, and we found complete agreement between the ap¬ 
proaches. With the exception of the pentagon integrals and M j we present the 
results for the master integrals as a Laurent expansion in the dimensional regulator up to 
terms of transcendental weight six for the soft and the hard regions and up to transcen¬ 
dental weight five for the collinear region. For M j)^ and M [4 we were only able to obtain 
the Laurent expansion up to terms of weight three and five respectively, which is however 
sufficient to compute the Higgs boson cross section up to finite terms 2 . Whenever we 
were able to do so, we also include results for the master integrals valid to all orders in e 
in terms of generalised hyper geometric functions, which can easily be expanded in e using 
the HypExp package [30]. Note that all integrals presented in this section satisfy recursion 
relations with respect to the space-time dimension [31]. We have checked that our results 
satisfy these dimensional recurrence relations. 


5.4 Analytic results for the hard region 


Mf 


1 

e(l-2 e) 

i + 2 + 4e + 8e 2 + 16e 3 + 32e 4 + 64e 5 + 0(e 6 ), 


(5.26) 


M -2 = - 4(3 4e jt (1 4e) 3 ^ 2 ( 1 , 1 , —e; 1 - e, 1 — 2e; 1) (5.27) 

18 96 1 1 

=-pH— q- + (-96 + I 2 C 2 ) + - (-64(2 + 6 OC 3 ) + 64C 2 - 320( s + I 86 C 4 

c c J e z e 

+ e (320( s - 992C 4 + 24C 2 Cs + 444C 5 ) + e 2 (992C 4 - 128C 2 Cs - 2368 (5 + 60C 3 2 
+ llll( 6 ) + 0 (e 3 ). 

5.5 Analytic results for the soft region 


M s_ r (4-4e)r(i-2e) 2 r(i + e) 

1 2 e 4 r (2 — 6e)r(l — e) 

= 4 - 4 + 4 (24 - I8C2) + - (112 + 24( 2 - 138( s ) + 672 - 144C 2 + 184( 3 
e * e 0 c z e 

+ (4032 - 672C 2 - 1104C 3 + 828( 4 + 828(2(3 - 4014( 5 ) e + (^24192 - 4032( 2 
- 4968(4 + 5352(5 - 1104( 2 (3 + 3174( 2 - ^^( 6 )e 2 + 0(e 3 ), 


(5.28) 

621( 4 

5152(3 


2 At least for the 37 first terms in the threshold expansion, see ref. [6]. 
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(5.29) 


M? = 


T(4 — 4e)T(l — 2e) 3 r(l + e) 
e 2 r(4 — 6e)T(l — e) 3 
1 11 61 / „ 55 A \ „ 305 A 484 869 

= Ms7 + T - 5C2 + ( 117 -T C2 - 446 ) e+ ( 687 -—— 

+ (4077 - 585C2 - + 22 OC 2 C 3 - 1332(4,) e 3 + (24327 - 3435< 2 

n 

e 4 + 0(e 5 ), 


~ 5148C3 - ^^( 4 + ^( 2(3 - 4884C5 + 968C 3 2 - 


16 


M .f = 


r(4 - 4e)T(2 - 3e) 

(1 — 2e) 2 er(4 — 6e)T(l — e) 

1 14 _ /400 


4 . /81920 


= j + y + (24 - 6 ( 2 )£ +{^~ 28 C 2 - 42C 3 ) 

+ (4576 - 8 OOC 2 - IOO 8 C 3 - 910C4 + 252C2C3 - 1302(g) - , v 

9219 \ 

- 4680C4 + 1176C2C3 - 6076C 5 + 882(1 - Ce)e 5 + 0 (e 6 ), 


2 /2320 

E +(— 

e 4 + 


(5.30) 

144( 2 - 196( 3 - 195( 4 )e 3 
- 4640(2 - 5600( 3 


Mf = 1 (-160 + 96(2) + - 704(2 + 672( 3 ) + ^(l536( 2 (5.31) 

- 4928(3 + 2436( 4 ) - 1024( 2 + 10752( 3 - 17864( 4 - 5760( 2 ( 3 + 16872( 5 
+ e (-7168(3 + 38976(4 + 42240( 2 ( 3 - 123728( 5 - 20160(f + 34710(e) + 0{e 2 ), 


Mi = - 


4r(4 - 4e)r(l - 3e) 


e(l + e)(l — 2e)T(3 — 6e)T(l — e) 
12 


3 7 ? 2 (1,1,1 — e; 2 — 3e, 2 + e; 1) 


(5.32) 


= -—( 2 - 8( 2 - 36(3 + (-112(2 - 24(3 + 33( 4 )e + (-672( 2 - 336( s + 22( 4 + 720( 2 ( 3 

/ 16881 \ 

- 450(5)e 2 + ( - 4032( 2 - 2016(3 + 308( 4 + 480( 2 ( 3 - 300( 5 + 1512( 3 2 + — — ( 6 je 3 

+ 0 ( 6 4 ). 



r(4-4e)r(l-2e) 3 r(l + e) 

2e 3 (l + 2e)r(l — 6e)T(l — e) 2 r(2 — e) 1 ' ’ 

= + 7 (~ 79 + 15 & + 161 “ 125 ^ 2 + 132(3 + ( - 319 + 395( 2 - 1100( 3 

+ — p -( 4 ) e + (641 - 805( 2 + 3476( 3 - ^^( 4 - 660( 2 ( 3 + 3996 ( 5 ) e 2 + ( - 1279 
+ 1595(2 - 7084( 3 + ^ Jp( 4 + 5500 ( 2(3 - 33300(5 - 2904( 3 2 + ^^( 6 ) e 3 + 0(e 4 ), 
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(5.34) 


M? = 


T(4 — 4e)r(l — 3e) 
e 3 (l - 2 e)r (2 - 6 e)r(l - e) 


3 F 2 ( l,-e,e;l- 2 e,l-e;l) 


6 


1 


= + -y + -(56 - 42C 2 ) + 336 - 28C 2 - 288( 3 + 2016 - 392( 2 - 192( 3 - 


2433 a> 


+ 


(l2096 - 2352( 2 - 2688( 3 - 811(4 + 1980( 2 ( 3 - 8262( 5 )e 2 + (72576 - 14112( 2 


204663 \ 

- 16128(3 - 11354(4 + 1320( 2 ( 3 - 5508( 5 + 6804(f--—( 6 )e 3 + 0(e 4 ) 


Mu = - 


2r(4-4e)r(l -2e) 3 r(l + e) 
e 3 T(3 — 6 e)r(l — e ) 3 


;F 2 (l,l-2e,e;2-3e,l-e;l) 


(5.35) 


6 16 1 

= —* - - 5 - + - (-112 + 36( 2 ) - 672 + 96( 2 + 306( s + (-4032 + 672( 2 + 816( 3 
+ 1410( 4 )e + (-24192 + 4032( 2 + 5712( 3 + 3760( 4 - 1842( 2 ( 3 + 8757( 5 )e 2 
+ (- 145152 + 34272(3 + 24192( 2 + 26320( 4 - 4912( 2 ( 3 + 23352( 5 - 7824(f 


+ 


57177 


6 3 + 0 (e 4 ), 


11 143 1 , _ 1 /4301 767 

M 9 — (-165 - 59( 2 ) + ^ + 


^(2 - 488( 3 ) + - ( - 5005 (5.36) 


+ 885( 2 + ^( 3 - ^( 4 ) + 15015 - 7320(3 + ^( 4 + 2576( 2 ( 3 

190808 „ 139365 _ 33488 


3 4 

- 14502(5 + e ( - 45045 + 26845( 2 - 

849273 \ . o. 

-(e) +0(e 2 ), 


(3 + 


-(4 - 
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5.6 Analytic results 
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6. Conclusions 

This paper discussed the expansion around threshold of the one-loop corrections to the 
production of a heavy colorless state in association with two partons. We introduced 
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techniques to compute the coefficients in the expansion, in principle to any desired order, 
and to express the result in terms of a small set of soft master integrals. These results 
are the missing pieces which went into the computation of the inclusive gluon-fusion Higgs 
production cross section as an expansion around threshold [6, 16, 19]. 

Our main tool to reduce the coefficients appearing in the threshold expansion to a 
linear combination of soft master integrals is reverse unitarity, which allows one to map 
phase-space integrals to cuts of Feynman integral. We perform this expansion separately 
for the phase-space measure and for the interference diagrams, and we observe that the 
phase-space integrals always reduce to integrals against the soft phase-space measure. The 
one-loop matrix elements are expanded in the soft limit using the strategy of regions, and 
the relevant regions are identified with the regions where the loop momentum is either 
hard, soft or collinear to one of the initial-state momenta. In each region, we combine the 
expanded interference diagrams with the corresponding phase-space measure and use IBP 
identities to reduce them to soft master integrals, independently of the region. 

The soft master integrals themselves are evaluated using two different approaches. 
The first method allows to derive a Mellin-Barnes representation for the soft integrals with 
poles at most at integer locations, provided that a similar Mellin-Barnes representation 
for the loop integration can be obtained. In our case most of the Mellin-Barnes integrals 
can be done in closed form to all orders in the dimensional regulator. In the remaining 
cases we were able to obtain a representation of the integral as a Laurent expansion in the 
dimensional regulator. The second method builds upon a specific factorisation of the phase 
space and exploits the knowledge of the scaling behavior of the integral with the external 
momenta to arrive at a one-fold parametric integral representation for the soft integral, 
which can be performed using modern integration techniques. 

The results for the soft master integrals presented in this paper are sufficient to obtain, 
at least, the first 37 terms in the threshold expansion of the N 3 LO gluon-fusion cross section, 
and conjecturally they provide the full set of boundary conditions to compute all the phase- 
space master integrals in general kinematics appearing in inclusive N 3 LO cross sections for 
the production of a heavy colorless state. We therefore anticipate that the results of this 
paper will have a substantial impact on future results for hadron collider cross sections at 
N 3 LO, e.g., Drell-Yan production or Higgs production via bottom-quark fusion. 
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A. Derivation of eq. (3.16) 

Here we wish to derive the tensor reduction of the following generic integral 


(P1,P2) 


d D k 

m D t 2 F(k,pi,p 2 ) 


(A.l) 


which was stated in eq. (3.16). A general tensor Ansatz would yield that the tensor integral 
j/u(pj ^p 2 ) can wr itten as a linear combination of {g^ v ,Pi,p%}■ But since 


k_L.pi = 0 = k±.p2, 


(A.2) 


we must also have that 


(Pi)ii k I lll " IJ ‘ n (PuP 2 ) = 0. for i = 1,2 and k = 1,2,.., n (A.3) 


and therefore I fll " tln (pi,P 2 ) can only depend on {g^, Pi, P%} through the transverse metric 


P 1 P 2 + P 2 P 1 

y± y 

P1-P2 


(A.4) 


Given that / /il "^ n (p 1 ,p 2 ) is fully symmetric under any permutations of the Lorentz indices 
pi ,.., p n , the tensor structure is fully determined to be that of eq. (3.16),i.e. 

I‘ xi ^ n (pi,p 2 ) = 7 ^g l i 1 ^ n I{n 1 pi,P2) (A.5) 

C ( n) 

where g^ " /i " is defined in eq. (3.19), C(n ) still has to be determined and 


I(n,PhP2) 


d D k ( k \) n / 2 
in 0 / 2 F(k,pi,P2) 


(A.6) 


To determine the coefficient C(n) we contract each of any n/2 of the n Lorentz indices on 
both sides of eq. (A.5) with any one of the remaining n/2 indices. Due to the symmetry it 
is irrelevant how this contraction is done and the result is 


C{n) = g^'^g^-gnu-ujin 


Using 

n 

Ul-.fln _ Ulfl2 U3--Vn , \ Ul £t2 M Ml • -l/k ■ -l/l ■ -M™ 

g± — g _l g± ■+■ /_j 

l^k=3 

we can explicitly carry out the contraction over p\ and P 2 to get 
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C(n) = (D ±g F "+ 2 Y, 

l>k=3 

= (D± + (n - 2 ))gf i “ lin g wlu ..g lln _ llMl 
= (Dj_ + (n — 2))C(n — 2) 


(A.7) 

(A.8) 


(A.9) 

(A.10) 

(A.ll) 
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(A.12) 


where D j_ = (g^) /J /t = D — 2 and in the second line we used the identity 

n 

” Ml-Mn _ 0 WMfc 

2 ^JL / ^_L ^_L 

Z>/c=l 

Given that C(2) = D_l it follows that 

n/2 

C{n) = f[(D ± + 2{i-l)) (A.13) 

2—1 

which completes the derivation of eq. (3.16). 
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